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Abstract 

We have carried out numerical simulation based on the equations of radiation 
magnetohydrodynamics to study the interaction of solar granules and small-scale 
magnetic fields in photospheric regions with various magnetic fluxes. Four sequences 
of 2D time-dependent models were calculated for photospheric regions with average 
vertical magnetic fluxes of 0, 10, 20, and 30 mT. The models exhibit no substantial 
variations in their temperature structure with varying average field strength, while 
the density and gas pressure profiles display gross changes. The solar granulation 
brightness field also varies substantially with magnetic flux. The contribution of the 
small-scale component to the intensity power spectrum increases with average field 
strength, whereas the large-scale component (of about a granule size) contributes 
less, the total rms intensity fluctuations being approximately the same. Thus the 
observed decrease in rms intensity fluctuations with growing average magnetic flux 
can be interpreted as smoothing of the small-scale component in the power spectrum 
by the modulation transfer function of the telescope. 

1 Introduction 

The interaction of solar granules and small-scale magnetic elements is of obvious interest, 
since about 90 percent of the total magnetic flux emerging from the solar photosphere is 
localized in these elements [21] . On the one hand, formation of small-scale elements is 
conditioned by the granulation scale of solar convection. On the other hand, the small- 
scale magnetic fields can have a profound impact on the solar granulation dynamics and 
structure. Thus theoretical models of small-scale concentrations should incorporate this 
scale of solar convection. Magnetic tubes are observed in the continuum as a network of 
bright points located in intergranular lanes [H] . The size of these points is beyond the 
resolution possibility of today's telescopes, so that the fine structure of these small-scale 
concentrations need be studied by combining theoretical methods [TSJ [19] with various 
observational data such as Stokes profiles of spectral lines [20], line asymmetries [21 E], 
observations with high spatial resolution [23] . 

By now a large number of semiempirical models of magnetic tubes have been de- 
veloped, they are able to account, at least qualitatively, for the structure of small-scale 
magnetic elements and to reproduce the observed Stokes profiles of solar spectral lines in 
magnetic regions. However, totally self-consistent models of small-scale magnetic elements 
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can be constructed only on the basis of the equations of radiation magnetohydrodynamics. 
Three-dimensional simulation of solar magnetoconvection was done by Nordlund [T5] and 
Nordlund and Stein [SI HZ], but the spatial resolution of their models is too small (the 
horizontal step of calculation grid is 100-200 km) to reproduce magnetic tubes of about 
100 km in size. Deinzer et al. [4,5], Knolker et al. |llj . Grossman-Doerth et al. [7] cal- 
culated detailed 2D models of small-scale magnetic tubes with a high spatial resolution. 
These models throw light on the dynamics and fine structure of small-scale magnetic con- 
centrations [TS], but they deal with single tubes only and ignore the interaction of tubes 
with each other. 

Here the object is to study the interaction between solar granulation and small-scale 
magnetic fields in photospheric regions with different magnetic fluxes. A detailed 3D 
simulation of this interaction requires an excessively great number of points in the cal- 
culation mesh and hence a prohibitively long machine time, so we used a complete set 
of 2D equations of radiation magnetohydrodynamics in the Cartesian coordinates for the 
solution of this problem. The basic idea was to construct solar granulation models with 
magnetic fields neglected and then to calculate a sequence of models with different den- 
sities of average magnetic flux with the same boundary and initial conditions (except for 
the initial magnetic field). In this way we could follow the influence of magnetic fields on 
granulation properties, brightness field, and fine structure of intergranular lanes, where 
magnetic tubes are formed. We note that we deal with magnetic fluxes on a plane (mag- 
netic slabs) in the 2D simulation, while these objects are magnetic flux tubes in actuality, 
they are simply called magnetic tubes. To avoid confusion, we use the term magnetic 
tubes in what follows. 



2 Radiation magnetohydrodynamics equations 

To simulate the interaction between solar granulation and small-scale magnetic fields, we 
used the complete set of 2D equations of radiation magnetohydrodynamics for stratified 
media in the Cartesian coordinates: 
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where i, j vary from 1 to 2, p and p are the medium density and gas pressure, are the 
velocity vector components, Bi are the components of the magnetic flux density vector, g is 
the vertical component of the gravitation vector, <$ i2 is the Kronecker delta, e = e + vf /2 
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is the sum of the internal and kinetic energy per unit mass, E = pe + B 2 /8tt is the 
total energy per unit volume, Qr is the divergence of radiant energy flux. The magnetic 
diffusion coefficient is a function of temperature, 77 = f]{T). The term Qj specifies the 
Joule dissipation: 

The gas pressure was found from the equation of state for a partially ionized gas: 



p = p(p,e). (6) 

The equation of state was calculated beforehand, and then, during calculations, the pres- 
sure p was found from tables. 

Turbulence effects were assumed to be simulated by numerical viscosity. 

3 Radiative transfer 

The quantity Qr, which specifies the radiative interaction between matter and radiation 
field in Eq. ([3]), was calculated from the relation 



Q R = 4tt J a v {B v - J v )dv, (7) 


where J v is the mean intensity, a v and B v are the absorption coefficient and the Planck 
function for the frequency v. To find Qr, we solved the radiative transfer equation for 
each time step: 

did did 

where f uxx and f uzz are the Eddington factors. The factors f uzz were determined from 
the 1-D solution of the radiative transfer equation along the OZ direction, while the 
factors f uxx were assumed to be 1/3 (this procedure is described in detail in [I]). Lateral 
boundary conditions were assumed to be periodic, and this suggested that the matter 
penetrated freely through the lateral boundary on one side and flew out with the same 
parameters on the opposite side. The upper and lower boundary conditions for Eq. (jHJ) 
were taken in accordance with Mihalas [13]. At the lower boundary, where the diffusion 
approximation is valid, we have 

dfuzzJu _ 1 dB v 

dz ~ 3^7' 1 j 

and at the upper boundary we have 



9 fuzz Jv 

dz 



a u h u J u . (10) 
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Here h v = H^/J^, H® being the vertical radiation flux at the upper boundary. The factor 
h v was found from the 1-D solution of the radiative transfer equation in the OZ direction. 
Equation flS]) was solved by the iteration method. For the simulation of upper photospheric 
layers to be correct, we allowed for the radiative transfer in spectral lines, as the energy is 
reradiated in this region predominantly in line frequencies. This was done with the OPDF 
table by Kurucz [12]. The whole frequency range was divided into four frequency groups 
depending on the total absorption coefficient (continuum + line opacities). The Planck 
function and the Rosseland mean absorption coefficient were tabulated for each frequency 
group. Thus radiative transfer equation (8) was first solved for four opacity groups and 
for each trime step, and then the integral quantity Qr was found with relation (J7|), each 
frequency range being taken with its weight. 



4 Calculation procedure 

To integrate the complete set of MHD equations (0Q)-(|ID> we used the conservative Total 
Variation Diminishing (TVD) scheme [8]. The principal concept of the scheme is to set 
limits on the general variation of numerical solution and prevent formation of nonphysical 
oscillations. The TVD scheme is second-order accurate in regions with a smooth solution, 
and it switches over to the first-order "against-the-flux" scheme near extrema. 
For the simple transfer equation 

where 

/ = au, (12) 
the TVD scheme can be written as in [6]: 



^^"-TI^-^)- (13) 



At 

where is the numerical flux 



Fi+i/2 = O-Hfi+i + fi) ~ 0.5|a m / 2 |[l - 4>(r i+l/2 )](u i+1 - m), (14) 

0(r i+1 / 2 ) is a switching function which depends on the ratio between the numerical gra- 
dients calculated in accordance with the flux direction: 

n+1/2 = , (15) 

Ui+i - Ui 

with er = sign(oti + i/2)- The switching function in this study was 

!2, when r > 2, 
r, when < r < 2, 
0, when r < 0. 

It switched TVD scheme (13)-(14) to a second-order scheme in regions with a smooth 
solution, 0(1) = 1, and to the first-order "against-the-flux" scheme near extrema, where 
0(r < 1) = 0. 

The second-order Adams-Bashforth scheme was used for the sampling of Eqs (l)-(4). 
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5 Boundary and initial conditions 



We assumed the lateral boundary conditions to be periodic. The boundary conditions 
at the upper boundary were taken open, but vertical velocities were scaled so that the 
mass flow through the upper boundary might be zero. Mean density and mean internal 
energy were specified by initial conditions, and their fluctuations were scaled so that the 
solution near the boundary might be smooth. A firm wall with a fixed internal energy 
was postulated at the bottom of the box simulated. We assumed that 

dB x , . 

ar = ° (16) 

for the magnetic field at the upper and the lower boundaries. 

The initial distribution of thermodynamic quantities was taken uniform at all hori- 
zontal levels, it was specified by the VAL80C model [21] for upper layers and by a model 
envelope calculated on the basis of the mixing-length theory [2] for deeper layers. A 
solenoidal velocity field in the form of a sum of harmonic functions, 

128 

u(x,z) = ^2Aisin(k xl x)ey.-p[-k zl (z max - z)], (17) 
i=i 



128 

w(x, z) = ^Aicos(k x ix)exp[-k z i(z max - z)], (18) 
i=i 

was introduced in the model to initiate convection in order to provide the greatest number 
of granulation motion scales in the box. 

The simulated region was represented by 256 x 128 cells of calculation mesh with the 
step Ax = Az = 15 km, it was rectangular in shape, 3840 km wide and 1920 km high. 
The upper boundary corresponded to the height h = 600 km above the tr = 1 level. 

The above algorithm was used to calculate a sequence of nonmagnetic models (with 
no magnetic field) for a simulation time of 16 min (which corresponds approximately to 
three passages of acoustic waves in the region). The matter in the box makes several 
revolutions in this time, and a statistically stable solution sets in. Then we introduced 
a vertical uniform magnetic field with different intensities into the model. Thus a non- 
magnetic model was employed as an initial condition for simulating magnetoconvection. 
The primary idea of the numerical experiment was to study changes in the solar granula- 
tion structure and properties in photospheric regions with different mean magnetic fluxes. 
With boundary and initial conditions specified, we have a single free parameter in the 
problem - the initial magnetic field intensity. In this numerical experiment we calculated 
four model sequences in the time interval from 16 to 25 min with initial mean magnetic 
field intensities of 0, 10, 20, 30 mT. 

6 The solar granulation dynamics and structure 
depending on magnetic flux 

Figure 1 shows the vertical and horizontal rms velocities in relation to the horizontally- 
averaged Rosseland optical depth for models with different magnetic fluxes. Our 2D 
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Figure 1: Vertical and horizontal rms velocities (V v and V/J as functions of the 
horizontally-averaged Rosseland optical depth for models with different mean magnetic 
fluxes: solid line) mT, dot-and-dash line) 10 mT, dashed line) 20 mT, circles) 30 mT. 

models being nonstationary, the velocities in Fig. 1 are averaged for each sequence of five 
models over 20-25 min intervals. As the initial magnetic field in our models is vertically 
directional, the mean field intensity has no pronounced effect on vertical velocities, while 
horizontal velocities are appreciably suppressed by the magnetic field in the upper layers, 
where the magnetic pressure exceeds or is comparable to the gas pressure. 

The simulation results do not reveal any substantial changes in the averaged structure 
of the models depending on magnetic flux, but the fluctuation profiles of thermodynamic 
quantities vary considerably with magnetic flux. Fluctuations in the total pressure in 
the upper layers diminish with growing magnetic flux due to the suppression of horizontal 
motions by the vertical magnetic field, i.e., the models tend to a quasi-equilibrium state in 
the horizontal direction with growing magnetic flux. At the same time relative fluctuations 
in the gas pressure and density increase with magnetic flux density owing to the growing 
magnetic pressure fluctuations. 

Figure 2a depicts the structure of the velocity field and temperature field for the 
models with a simulation time of 25 min and different mean magnetic fluxes. Temperature 
inhomogeneities of about granule size become more smoothed with increasing mean field 
intensity, while the fluctuations in intergranular lanes grow. Temperature fluctuations in 
the optically thin upper layers are governed primarily by radiative transfer, that is why 
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Figure 2: Results of simulation: a) isotherms (lines) and velocity field (arrows); thick 
line shows a 6400 K isotherm, which corresponds approximately to the level Tr — 1 for 
homogeneous models; b) isopycnics (thick lines) and lines of magnetic field force (thin 
lines). The simulation time is 25 min. Mean vertical magnetic field intensity is 0, 10, 20, 
and 30 mT (from top to bottom). 
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the temperature fluctuations at small optical depths are close in magnitude in all models. 
The area occupied by downflows increases as well, and regions where the velocities are 
suppressed by magnetic fields appear in intergranular lanes in the models with a mean 
magnetic flux of 20 mT or higher. 

Figure 2b shows isopycnics together with lines of magnetic force. The solar convec- 
tion is known to gather the lines of magnetic force together and form vertical bundles 
in intergranular lanes. The vertical structure of magnetic tubes is lost near the base of 
the region simulated, since the convective cells are closed in this region because of closed 
boundary conditions, and the velocities become predominantly horizontal, thus leading 
to the destruction of vertical structure of magnetic tubes. The magnetic tubes expand 
in the upper photosphere due to decreasing pressure with growing height, they become 
thicker, and the area occupied by them grows with the mean field intensity. While the 
model temperature structure is not markedly affected by variations in the mean magnetic 
flux, the pressure distribution changes drastically. Abrupt density drops appear in inter- 
granular lanes, the widths of these drops grows with the width of magnetic tubes, and 
magnetic tubes occupy a considerable part of the area in models with a mean intensity 
of 30 mT. Density profiles above granules become more sinuous. 

7 Fine structure of magnetic tubes 

Figure 3 illustrates the structure of a magnetic tube on an enlarged scale. Isopycnics and 
lines of magnetic force are drawn together with the velocity field for geometric area 2 
shown by a square in Fig. 2b. Drops in density profiles grow with magnetic flux. Motions 
in intergranular lanes are suppressed by the magnetic field. Downward streams flow 
around the magnetic tube, and so the area occupied by them increases. An interesting 
effect can be seen in the area with a 30 mT mean intensity. The magnetic tube becomes 
thicker and breaks up into two tubes. Motions are found to be suppressed in the space 
between two tubes, while there are downflows inside the tubes. This is likely to be a 
structure intermediate between a facular point and a pore; a similar configuration was 
simulated by Knolker et al. [10]. Unlike in the magnetic tube in region 2, motions in the 
tube lying in geometric region 1 (Fig. 2b) are completely suppressed by the magnetic field 
at a mean intensity greater than 10 mT. This tube resembles the classical magnetic tube 
model. 

The surface of equal gas and total (gas + magnetic) pressures for the geometric regions 
shown by squares 1 and 2 in Fig. 2b is a tube-like structure of small-scale magnetic ele- 
ments. The width of drops in the gas pressure profile grows with magnetic field intensity. 
Two drops in the gas pressure profile can be seen in the models with a mean intensity 
of 30 mT, where a tube breaks up into two magnetic tubes (Fig. 3d), but the surface of 
equal total pressure does not vary with magnetic flux. The total pressure remains ap- 
proximately constant at all horizontal levels, i.e., the models tend to a quasi-equilibrium 
state in the horizontal direction while the mean magnetic flux increases. 

8 Brightness variations depending on magnetic flux 

Continuum intensity profiles calculated for models with different magnetic fluxes for 
A = 500 nm and for a simulation time of 25 min are plotted in Fig. 4. The model 
with no magnetic field displays bright granules 1000-2000 km in size. The granule bright- 
ness amplitude decreases in the models with a magnetic field, but brightness peaks with 
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Figure 3: Isopycnics (thick lines), lines of magnetic field force (thin lines), and velocity 
field for the geometric region shown by square 2 in Fig. 2b on an enlarged scale. Mean 
vertical magnetic field intensity is 0(a), 10(b), 20(c), and 30(d) mT. 

intensities up to 1.6 appear in intergranular lanes. Deinzer et al. |5j point out that these 
intensity peaks are surrounded by dark rings formed by local peaks in the density profiles 
near magnetic tube boundaries (Fig. 3). The area occupied by bright points increases 
with magnetic flux, and the amplitude of large-scale inhomogeneities (about a granule 
size) decreases. Overall fluctuations of the intensity averaged over five models change 
only slightly with varying magnetic field intensity, they are about 14% for A = 500 nm. 

To demonstrate the intensity field "smearing" due to a limited spatial resolution, we 
convoluted the intensity profiles calculated with 2D models with a modulation transfer 
function (MTF) for an ideal telescope of diameter D = 30 cm. This function was rep- 
resented by a Gaussian profile with the half- width X/D. Figure 4 shows the intensity 
profiles smoothed in this manner. The intensity contrast decreases considerably, and 
bright points are hardly discernible among granules. Besides, the horizontal dimension of 
smoothed bright peaks in models with a large mean magnetic flux becomes comparable 
to the size of small granules. Therefore, it is difficult to distinguish the bright points 
associated with small-scale magnetic structures from granules in facular plages with large 
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Figure 4: Radiation intensity in the continuum for models with a 25-min simulation 
time and different mean magnetic fluxes (solid line). Dashed line shows intensity profiles 
obtained after the convolution with the modulation transfer function of the telescope. 
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Figure 5: Root-mean-square fluc- 
tuations in the continuum in- 
tensity calculated for 2-D mod- 
els (solid line with circles) and 
smoothed by the modulation 
transfer function of the tele- 
scope (dashed line with circles); 
original observations (solid line 
with squares) and observations 
with 5-min oscillations filtered off 
(dashed line with squares) from 
the data of |22|. 
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magnetic filling factors. 

The intensity power spectrum averaged over five models with different magnetic fluxes 
suggests that the small-scale component grows with field intensity while the large-scale 
component becomes weaker. When observed with a limited spatial resolution, small-scale 
bright points are smoothed by the telescope MTF, and the observed overall fluctuations 
of intensity diminish with growing mean magnetic flux. 

Figure 5 shows the intensity fluctuations calculated from 2D models and observed 
with the SOUP space telescope with a 30 cm aperture [22] ■ As there is no fundamental 
difference between the positive and negative orientations of the initial magnetic field in our 
models, the intensity fluctuations were assumed to be symmetric about the zero magnetic 
field. Obviously the intensity fluctuations calculated from 2D models are approximately at 
the same level for the field intensity range from to 30 mT, but the fluctuations smoothed 
by the telescope MTF diminish with growing mean field intensity. The observations in [22] 
display the same tendency. Thus, the observed decrease of intensity fluctuations in the 
continuum with increasing mean magnetic flux can be accounted for, at least qualitatively, 
within the framework of our 2D models. 

9 The Wilson depression 

The Wilson depression (Fig. 6) is about 150-200 km in the regions of magnetic flux 
concentration. Bright peaks arise in the continuum, as the horizontal temperature struc- 
ture varies more smoothly than the height of the layer where the optical depth is unity. 
Besides, such depressions of the level tr — 1 are outlets for the radiant energy from hot 
lower layers, which results in additional heating of the upper parts of magnetic tubes [TTj . 
The relationship between magnetic flux density and distance (Fig. 6) suggests that there 
is a clear correlation between magnetic flux in tubes and bright intensity peaks. There is 
also a clear distinction between small-scale magnetic elements and nonmagnetic regions 
with a nearly zero magnetic field intensity. The width of magnetic tubes grows with 
field intensity, but the maximum magnetic flux density in small-scale magnetic structures 
remains nearly constant and is about 200 mT. 

So, we constructed four sequences of model solar photospheres with different magnetic 
fluxes. When a magnetic field is introduced in the models, three main free parameters 
appear: initial magnetic field configuration, magnetic diffusion, and initial field intensity. 
The initial magnetic field orientation, which determines the predominant direction of 
magnetic flux, was fixed in our simulation (the initial magnetic field was assumed to 
be vertical). The magnetic diffusion coefficient specifies the maximum concentration of 
magnetic flux density and, consequently, other parameters that depend on these quantities 
such as drops in pressure and density, the Wilson depression in the magnetic tube. The 
magnetic diffusion coefficient was a function of temperature in the simulation and was 
defined by one and the same law for all models. In the calculations we chose the smallest 
magnetic diffusion coefficient that ensured a stable numerical solution, it was about rj = 
10 16 cm 2 s _1 at T = 10000 K. Thus, only the initial field intensity was varied in our 
models. The maximum magnetic flux density being limited by the magnetic diffusion 
coefficient, a growth of the mean magnetic flux resulted in widening of magnetic tube, 
i.e., in an increase of magnetic filling factor. The filling factor for MHD models should be 
estimated, however, from 3D calculations, since the topologies of 2D and 3D convective 
motions are fundamentally different. 
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Figure 6: Geometric height for the optical depth tr = 1 (a) and vertical magnetic flux 
density at the level tr = 1 (b) v. horizontal distance for models with different mean 
magnetic fluxes. 

10 Conclusion 

The major results of our study follows. 

1. A vertical magnetic field with a mean intensity of 10-30 mT has no appreciable 
effect on vertical velocities, while horizontal velocities are suppressed to a considerable 
degree by the field. 

2. Fluctuations in the total pressure diminish with growing mean field intensity owing 
to the suppression of horizontal velocities by a vertical magnetic field, i.e., the models tend 
to a quasi-equilibrium state in the horizontal direction when the magnetic flux increases. 

3. Fluctuations in the relative gas pressure and density grow with field intensity, while 
the temperature fluctuations grow in subphotospheric layers only and become approxi- 
mately the same in all models at the upper levels. 

4. Gas motions are suppressed inside magnetic tubes in the models with a mean 
intensity of 20 mT and higher. 

5. A peculiar configuration appears in the models with a mean field intensity of 30 mT: 
a thick tube breaks up into two tubes, with motions suppressed in the region between the 
tubes rather than inside the tubes. This is likely to be a structure intermediate between 
a facular point and a pore. 

6. A magnetic field affects only slightly the overall fluctuations of the continuum 
intensity, but changes substantially the intensity power spectrum. The contribution of the 
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small-scale component to the intensity power spectrum grows with magnetic flux, while 
the contribution from the large-scale component (of about a granule size) diminishes, so 
that observations with a limited spatial resolution show a decrease of intensity fluctuations 
with growing magnetic flux, as small-scale bright spots are smoothed by the modulation 
transfer function of the telescope. 

The approach elaborated here has some limitations and disadvantages. The two- 
dimensional representation of granulation motions is the crudest assumption in our treat- 
ment. Magnetic tubes in our models are formed between two granules, while small-scale 
elements in the solar photosphere are observed in the lanes formed by at least three or 
more neighboring granules. We dealt with the granular scale of solar convection only, and 
the initial magnetic field was artificially introduced in the models. Formation of initial 
field need be studied before the formation of facular plages and network areas can be 
understood. This problem can be studied by simulating larger solar convection scales 
such as meso- and supergranulation. 

In spite of the shortcomings of our method, the 2D models calculated by us can be 
used for the interpretation of observational data for solar photospheric regions with various 
mean magnetic fluxes. The models will be applied in our subsequent studies to calculate 
profiles of actual solar photospheric lines. 

Acknowledgments. We wish to thank A. S. Gadun for useful discussion of the 
results. 
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